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ABSTRACT 

Due to the chaotic nature of the Solar System, the question of its long-term stability can 
only be answered in a statistical sense, for instance, based on numerical ensemble integrations of 
nearby orbits. Destabilization of the inner planets, leading to close encounters and/or collisions 
can be initiated through a large increase in Mercury’s eccentricity, with a currently assumed 
likelihood of ^1%. However, little is known at present about the robustness of this number. Here 
I report ensemble integrations of the full equations of motion of the eight planets and Pluto over 
5 Gyr, including contributions from general relativity. The results show that different numerical 
algorithms lead to statistically different results for the evolution of Mercury’s eccentricity (ex). 
For instance, starting at present initial conditions (ex — 0-21), Mercury’s maximum eccentricity 
achieved over 5 Gyr is on average significantly higher in symplectic ensemble integrations using 
heliocentric than Jacobi coordinates and stricter error control. In contrast, starting at a possible 
future configuration (ex — 0.53), Mercury’s maximum eccentricity achieved over the subsequent 
500 Myr is on average significantly lower using heliocentric than Jacobi coordinates. For example, 
the probability for ex to increase beyond 0.53 over 500 Myr is >90% (Jacobi) vs. only 40-55% 
(heliocentric). This poses a dilemma as the physical evolution of the real system — and its 
probabilistic behavior — cannot depend on the coordinate system or numerical algorithm chosen 
to describe it. Some tests of the numerical algorithms suggest that symplectic integrators using 
heliocentric coordinates underestimate the odds for destabilization of Mercury’s orbit at high 
initial ex- 

Subject headings: celestial mechanics — methods: numerical — methods: statistical — planets and 
satellites: dynamical evolution and stability 


1. Introduction 

The question whether the Solar System is dy¬ 
namically stable over long periods of time has 
received considerable attention for centuries, in¬ 
cluding contributions from Newton, Lagrange, 
Laplace, Poincare, Kolmogorov, Arnold, Moser 
etc. (e.g. Laskar 2013). Research in this field 
has recently experienced a renaissance due to ad¬ 
vances in computational power and numerical al¬ 


gorithms, permitting integration of the full equa¬ 
tions of motion approaching the Solar System’s 
lifetime (± ^5 Gyr) (Wisdom & Holman 1991; 
Quinn et al. 1991; Sussman & Wisdom 1992; Saha 
& Tremaine 1992; Murray & Holman 1999; Ito 
& Tanikawa 2002; Varadi et al. 2003; Batygin & 
Laughlin 2008; Laskar & Gastineau 2009). More¬ 
over, parallel computing now allows tackling the 
problem of stability with statistical means through 
simultaneous integration of multiple nearby orbits 
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(Laskar & Gastineau 2009). A statistical approach 
is necessary because of the chaotic behavior of the 
system, i.e. the sensitivity of orbital solutions to 
initial conditions (Laskar 1989; Sussman & Wis¬ 
dom 1992; Murray & Holman 1999; Richter 2001; 
Varadi et al. 2003; Batygin & Laughlin 2008). 
Small differences in initial conditions grow expo¬ 
nentially, with a time constant (Lyapunov time) 
for the inner planets of only ~5 Myr (Laskar 1989; 
Varadi et al. 2003; Batygin & Laughlin 2008). For 
example, a difference in initial coordinates of 1 mm 
grows to ~1 AU (= 1.496x10^^ m) after 163 Myr. 
Thus, the chaotic nature of the planetary orbits 
makes it fundamentally impossible to predict their 
evolution accurately beyond ~100 Myr (Laskar 
1989). Hence the stability problem can only be 
answered in probabilistic terms, e.g. by studying 
the behavior of a large number of physically possi¬ 
ble solutions. The quest for a single deterministic 
solution, which conclusively describes the Solar 
System’s evolution for all times (in the spirit of 
Laplace’s demon, Laplace 1951), must be regarded 
as quixotic. 

An ensemble integration of 2,501 nearby orbits 
over 5 Gyr has recently been reported based on 
the full equations of motion and including con¬ 
tributions from general relativity and the Moon 
(Laskar & Gastineau 2009). In those simula¬ 
tions, two adjacent orbits differed initially by 
only 0.38 mm in Mercury’s semi-major axis (a^vi). 
The largest overall offset in initial gm was hence 
2, 500 X 0.38 mm = 0.95 m, well within the un¬ 
certainty of our current knowledge of the Solar 
System. In about 1 % of the simulations, a large 
increase in Mercury’s eccentricity was reported, 
which can lead to destabilization of the inner plan¬ 
ets, including close encounters and/or collisions. 
Note that the obtained probability distribution 
of Mercury’s eccentricity was similar to results 
obtained earlier with averaged equations (Laskar 
2008). Most surprisingly, further simulations in¬ 
cluded the possibility of a collision between Earth 
and Venus via transfer of angular momentum from 
the giant planets to the terrestrial planets (Laskar 
& Gastineau 2009). Given the chaotic behavior 
of the system, such an outcome might be within 
the range of possibilities. However, little is cur¬ 
rently known about the potential dependence of 
the simulated trajectories and statistical results 
on the numerical integrator, step size, integration 


coordinates, etc. 

Here I report ensemble integrations of Solar 
System orbits over 5 Gyr, including contributions 
from general relativity. The simulations reveal a 
strong influence of integration coordinates on sta¬ 
tistical results, i.e. on the predicted odds for desta¬ 
bilization of Mercury’s orbit. Furthermore, I dis¬ 
cuss resonances and Lyapunov times, and perform 
several tests aiming at resolving the dilemma of 
the dependency of statistical results on integration 
coordinates. 

2. Methods 

The full equations of motion of the eight plan¬ 
ets and Pluto were integrated over 5 Gyr into 
the future using the numerical integrator pack¬ 
ages mercuryS and HNBody and various integra¬ 
tion options (Chambers 1999; Rauch & Hamil¬ 
ton 2002). Unless stated otherwise, symplectic 
integrators were used. Symplectic integrators ex¬ 
actly describe the time evolution of a (modified) 
Hamiltonian system that is very close to the true 
Hamiltonian (cf. e.g. Wisdom & Holman 1991; 
Yoshida 1993; Chambers 1999) and hence do not 
suffer from long-term buildup in energy error. The 
maximum energy variation along a given orbit 
in symplectic integrations (see below) provides a 
measure of the difference between the modified 
and the true Hamiltonian. Relativistic corrections 
are critical (Laskar & Gastineau 2009) and were 
available in HNBody but not in mercuryG. Post- 
Newtonian corrections for symplectic integration 
(Mikkola 1998; Soffel 1989) were therefore imple¬ 
mented before using mercuryG (see Appendix). 
Thus, all simulations presented here include con¬ 
tributions from general relativity, unless explicitly 
stated. To allow comparison with previous stud¬ 
ies, several higher-order effects were not included 
here. This means ignoring potential effects from 
asteroids (Ito & Tanikawa 2002; Batygin & Laugh¬ 
lin 2008), perturbations from passing stars, and 
solar mass loss (Ito & Tanikawa 2002; Batygin & 
Laughlin 2008; Laskar & Gastineau 2009). Fur¬ 
thermore, the Earth-Moon system was considered 
a single mass point, located at the Earth-Moon 
barycenter. 

Per ensemble integration, 40 orbital solutions 
were computed with each package, starting from 
the same set of initial conditions, where Mercury’s 
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initial radial distance was offset by 7 cm between 
every two adjacent orbits (largest overall offset: 
39 X 7 cm = 2.73 m). Initial conditions for all bod¬ 
ies in the 5-Gyr runs (before offsetting Mercury) 
were generated from DE431 (naif . jpl .nasa. 
gov/pub/naif/generic_kernels/spk/planets) 
at JD2451544.5 (01 Jan 2000) using the SPICE 
toolkit for Matlab (naif.jpl.nasa.gov/naif/ 
toolkit.html) (Table 2). 

For the 5-Gyr simulations with mercuryS, the 
hybrid integrator with democratic-heliocentric co¬ 
ordinates (Chambers 1999; Duncan et al. 1998) 
following Batygin & Laughlin (2008) and a 6-day 
initial timestep (At) for the 2nd order symplec- 
tic algorithm was used. In case Mercury’s eccen¬ 
tricity (cm) increased above certain threshold val¬ 
ues during the 5-Gyr simulations with mercuryG, 
At was reduced but held constant after that un¬ 
til the next threshold was crossed (if applica¬ 
ble). For example, starting at At = 6 days and 
Cm = 0.21, At was reduced to 5 days when bm 
exceeded 0.34. mercuryG’s source code is avail¬ 
able (Chambers 1999) and was manipulated ac¬ 
cordingly. The threshold values and correspond¬ 
ing timesteps (Table 1) were chosen to maintain 
a roughly similar maximum relative energy error 
(max|A£'/£'| = max|(E(t) — £'o)/Eo|) throughout 
the simulation. HNBody’s source code is not avail¬ 
able (Rauch & Hamilton 2002) and a constant 4- 
day timestep was used throughout the 5-Gyr sim¬ 
ulations, employing a 2nd order symplectic inte¬ 
grator with corrector and Jacobi coordinates (see 
below and Wisdom & Holman 1991). 


Table 1: Threshold cm values and timesteps for 
the 5-Gyr simulations with mercuryG. 


Cm > 

0.34 

0.38 

0.44 

0.49 

0.57 

0.62 

At (d) 

5 

4 

3 

2 

1.5 

1 


The 40 orbital solutions were simultaneously 
integrated on a 64-bit Linux cluster (one 5-Gyr 
job per core on 10 Intel 17-3770 3.40 GHz nodes 
with 4 cores each). Typical integration times for 
the 5-Gyr runs with mercuryG were 12 — 23 days 
wall-clock time depending on step-size reduction 
{At < 6 d, see above) and 12 days wall-clock time 
with HNBody {At = 4 d). Integrations with HNBody 
were performed using the pre-compiled binary ver¬ 


sion for Linux 64-bit Intel/AMD64. Executa¬ 
bles of mercuryG were compiled from source code 
on 64-bit Linux platforms using gf ortran 4. G . 3. 
When compared to test runs with mercuryG on a 
true 32-bit machine, a small long-term decrease in 
angular momentum was noticeable in the output 
of the 64-bit mercuryG executable at small step 
size. The issue disappeared when using one of the 
following options on 64-bit platforms. (1) Com¬ 
pile with flag -m32 (generates code for a 32-bit 
environment but is slower). (2) In the mercuryG 
subroutine drift_kepmd() replace the original ap¬ 
proximations for sin(a;) and cos(a;) by actual sin(a;) 
and cos(a;) fnnctions. 

2.1. Democratic-heliocentric- and Jacobi 
coordinates 

The effect of different integration coordinates 
on the results of the symplectic integrations will 
be discussed at length below. Hence the defini¬ 
tion of the coordinates is given here. Democratic- 
heliocentric coordinates simply consist of heliocen¬ 
tric (bodycentric) positions and barycentric mo¬ 
menta (Duncan et al. 1998; Chambers 1999), see 
Section 6.1. 

Jacobi coordinates have a slightly more com¬ 
plex structure but are useful in symplectic inte¬ 
grations of the n-body problem because they al¬ 
low writing the kinetic energy as a diagonal sum 
of squares of the momenta (Wisdom & Holman 
1991). In Jacobi coordinates, the first coordinate, 
Qq, may be taken as the position of the center 
of mass. The first relative coordinate is the posi¬ 
tion of the first planet relative to the central body. 
The second relative coordinate is the position of 
the second planet relative to the center of mass of 
the central body and the first planet, and so forth. 
Thus, the ith relative coordinate is the position of 
the ith planet relative to the center of mass of the 
central body and the planets with lower indices: 


Qi — ^i—1 

where 0 < z < n and: 

(1) 

i 

Xi = ^^rrijXj ; 

II 


j=o j=o 


is the center of mass of bodies with indices up to 
i. The momenta are: 

p^ = m'y, , (3) 
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where Vi is the time derivative of Qi and m' = 
Mi-irrii/Mi for 0 < i < n and rriQ = M„_i = M, 
is the sum of all masses. Interestingly, because of 
the Sun’s dominant mass, the difference in the ac¬ 
tual values of heliocentric and Jacobi coordinates 
is small in the Solar System (e.g. zero for Mer¬ 
cury and of order 10“^ and 10“® AU for Venus 
and Earth). Yet the difference in the results of 
symplectic integrations due to these different sets 
of coordinates is significant (see below). Note 
that the comment above is meant for illustrative 
purposes only. All initial input coordinates for 
mercuryS and HNBody were specified in heliocen¬ 
tric positions and velocities (e.g. Table 2). Trans¬ 
formation to the desired set of integration coor¬ 
dinates is performed internally. In the following, 
’’same initial conditions” refers to truly identical 
initial conditions, regardless of the integration co¬ 
ordinates used. 

3. Results: 5-Gyr and 500-Myr simula¬ 
tions 

At first glance, the 5-Gyr simulations with 
mercuryS and HNBody appeared to yield similar 
results, except for different energy and momen¬ 
tum errors (Fig. 1). In none of the 80 simu¬ 
lations a large increase in cm or a destabiliza¬ 
tion of the inner solar system occurred, which is 
also not expected if the corresponding odds are 
~1% (Laskar & Gastineau 2009). However, at 
second glance differences between the mercuryB 
and HNBody runs became apparent. For runs with 
identical initial conditions, the maximum differ¬ 
ence in Mercury’s eccentricity (Abm ) between the 
mercuryS and HNBody simulations typically grew 
from ^10“^ to 0.01 over only ^10 Myr. For com¬ 
parison, when HNBody was run with bodycentric 
coordinates (all else equal), Ae^vi grew to 0.01 
over ^60 Myr (see below). Furthermore, Ae^vi 
between two simulations with the same package 
and options but different initial conditions grew 
to 0.01 over ^75 Myr and ^130 Myr, respec¬ 
tively, for the above mercuryG and HNBody setup 
(Fig. 5). While the divergence of trajectories after 
^60 Myr may be attributed to the chaotic behav¬ 
ior of the physical system (Laskar 1989; Varadi 
et al. 2003; Batygin & Laughlin 2008), the rapid 
rise in Aex over only ~10 Myr between the 5- 
Gyr runs with democratic-heliocentric (mercuryG) 
and Jacobi coordinates (HNBody) hints to a po¬ 





Time (Gyr) 


Fig. 1.- Results from 5-Gyr runs with mercury6 and 

HNBody. Blue curves (a) —(d): mercury6, run ^2 (solution 
5^®), At < 6 days, democratic-heliocentric coordi¬ 

nates. Green curves (a) —(d): HNBody run #5 
At — 4 days, Jacobi coordinates, (a) Relative energy error, 
\AE/E\ — \{E{t) — Eo)/Eq\. (b) Relative angular momentum 
error, \AL/L\ — \{L{t) — Lq)/Lq\. (c) Change in Mercury’s 
semi-major axis, Aa = a(t) — uq. (d) Mercury’s eccentricity 
(s) Mercury’s maximum eccentricity (per 1 Myr bin) 
from all 5-Gyr runs {N — 80) with mercuryB and HNBody 
(AT = 40 each). Note time step reduction from 6 to 5 days at 
'^2.24 Gyr in (see (a) and (c)). 


tential numerical origin. A rapid Ae^vi rise over 
^10 Myr was also observed between two HNBody 
runs with identical initial conditions and timestep 
but bodycentric vs. Jacobi coordinates. In addi¬ 
tion, orbital solutions obtained with bodycentric 
and Jacobi coordinates yield different Lyapunov 
times for the inner planets (see below). 

The maximum in the relative energy error in 
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the 5-Gyr mercuryS runs was typically 10^ times 
larger than in HNBody; the corresponding maxi¬ 
mum angular momentum error about 10^ times 
larger (Fig. 1). One might hence be tempted to 
argue that any differences between the present 5- 
Gyr simulations are due to larger errors in the 
mercuryG simulations (less reliable), rather than 
heliocentric vs. Jacobi coordinates. However, this 
argument is not supported by the results of ad¬ 
ditional ensemble simulations (see below). The 
small offsets in Mercury’s initial radial distance 
randomized the initial conditions and led to com¬ 
plete divergence of eccentricities after ^100 Myr 
(Fig. 1). No patterns in the evolution of Mer¬ 
cury’s eccentricity were observed between adjacent 
orbits (or sets of orbits); neither appeared pairs of 
mercuryG and HNBody runs starting from identical 
initial conditions lead to correlations in the behav¬ 
ior of eM over 5 Gyr. 

However, statistically the 5-Gyr ensemble sim¬ 
ulations with mercuryG (heliocentric coordinates) 
and HNBody (Jacobi coordinates) led to signifi¬ 
cantly different distributions for the maximum ec¬ 
centricities of Mercury (e(^’‘) achieved over 5 Gyr 
(iV = 40 each. Fig. 2). More precisely, the null hy¬ 
pothesis that the two samples from the 5-Gyr 
mercuryG and HNBody ensemble simulations are 
random samples from normal distributions with 
equal means can be rejected {p < 0.017, two-tailed 
Student’s t-test). Furthermore, the probability for 
Cm to grow, say, beyond 0.4 is more than double 
for the 5-Gyr mercuryG setup (11/40 = 28%) than 
for the HNBody setup (5/40 = 13%) (Fig. 2). 

This statistical discrepancy (obtained at rel¬ 
atively low initial bm — 0.21) raises questions 
about algorithm performance at high cm- For 
instance, one might expect that cm also grows 
larger on average at initially high bm when using 
heliocentric vs. Jacobi coordinates. The issue is 
critical for the potential destabilization of the in¬ 
ner planets as a result of a large increase in Mer¬ 
cury’s eccentricity. Thus, I conducted additional 
ensemble integrations over 500 Myr starting at the 
time/conditions of the run with the highest Bm — 
0.53 (solution achieved during the 5-Gyr sim¬ 
ulations (Figs. 2 and 3). The timestep was reduced 
to 2 days in both the mercuryG and HNBody setup. 
Surprisingly, in this case, the mercuryG simula¬ 
tions (heliocentric coordinates) gave a significantly 
smaller mean b’J^^ value than the HNBody simula- 





Fig. 2. - Mercury’s maximum eccentricity achieved 

during 5-Gyr and 500-Myr runs with mercury6 (blue bars) 
and HNBody (green bars). Bod: bodycentric, Jac: Jacobi 
coordinates. (a) Mean = 0.377. ej\yi grows beyond 

0.4 in 11 out of 40 solutions (i?o .4 — 11/40 — 28%). (b) 
Mean = 0.349; Rq.a ^ 5/40 ^ 13%. (c) Mean 

= 0.535; ej^ grows beyond start value of 0.53 

in 17 out of 40 solutions (i^■|■ — 17/40 — 43%). (d) 

Mean ^ 0.557; R-^ ^ 37/40 ^ 93%. (e) Mean 

^ 0.538; R^ - 21/40 - 53%. (f) Mean - 0.536; 

= 22/40 = 55%. The null hypothesis that two pairwise 
samples as shown in panels a. — a, c, e, f and 0 — b, d 
(^N — 40 each) are random samples from normal distributions 
with equal means can be rejected at the following signifi¬ 
cance levels Pa 0 (two-tailed Student’s t-test): pab < 0.017, 
Pcd < 2.2x10”'’, Ped < 0.0003. Pfd < l.SxlO”". 


tions (Jacobi coordinates) {p < 2.2x10“®). Also, 
the probability for bm to increase above the start 
value of 0.53 is less than half for the mercuryG 
setup (17/40 = 43%) than for the HNBody setup 
(37/40 = 93%) (Fig. 2). 

This result appears unrelated to the integrator 
package but related to the choice of the integration 
coordinates. Using HNBody with bodycentric coor¬ 
dinates for the 500-Myr runs gives a mean Bj^^^ 
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Fig. 3.— Example of Mercury’s eccentricity (e^) during a 
500-Myr run with HNBody (At = 2 days, Jacobi coordinates). 
The pattern is typical during increases associated with the 
(pi—ps) resonance (see text). 

value similar to the heliocentric mercuryS setup 
and significantly lower than the HNBody setup with 
Jacobi coordinates {p < 0.0003, Fig. 2). I also 
tested whether the statistically different at 

Mercury’s higher initial eccentricity is related to 
typically larger errors associated with bodycentric 
vs. Jacobi coordinates at the same step size. I 
repeated the 500-Myr runs with the bodycentric 
HNBody setup but an eight-fold smaller timestep 
(0.25 d, ^100-fold smaller \AE/E\). However, the 
basic result remained the same. At high initial 
Cm , the setup using bodycentric coordinates leads 
to significantly smaller mean than the setup 
using Jacobi coordinates (p < 1.3x10“®, Fig. 2). 
This suggests that the statistical discrepancies are 
due to integration coordinates, not errors in en¬ 
ergy or angular momentum (see above). 

This poses a quandary because the evolution 
of the physical Solar System, including its proba¬ 
bilistic behavior, cannot depend on our choice of 
coordinate system or numerical algorithm. Fun¬ 
damentally, there is no reason to prefer one set of 
coordinates over the other. However, symplectic 
schemes differ in their implementation of differ¬ 
ent integration coordinates (Section 6.1). While 
in symplectic algorithms the size of the pertur¬ 
bation term (planet interactions) may be numer¬ 
ically larger in heliocentric than in Jacobi coor¬ 
dinates, little performance difference has been re¬ 
ported over 100,000 steps (Farres et al. 2013). Fur¬ 
thermore, heliocentric and Jacobi coordinates ap¬ 
pear to have opposing effects on at low and 
high initial e^, pointing to a more complex origin. 
Importantly, differences per timestep due to inte¬ 
gration coordinates, which may ultimately cause 
differences in eccentricity, are minuscule. 


416 417 418 419 420 


For example, in HNBody (At = 6 days, body vs. 
Jacobi) it takes ^200 kyr, or 12 million steps for 
Acm to grow from 10 “® to 10 “^ (^^1 billion steps 
for Jupiter). Statistical differences between 
in the 500-Myr runs only become significant (p < 
0.05) after ~140 Myr, or 25 billion steps (HNBody, 
At = 2 days, body vs. Jacobi). 

4. Fourier analysis: 51—55 Resonance 

Fourier analysis of Mercury’s longitude of peri¬ 
helion showed that Mercury’s eccentricity increase 
was generally associated with a shift in eigenfre- 
quency ( 51 ) towards Jupiter’s forcing frequency 
( 55 ). For most 500-Myr runs, a correlation be¬ 
tween and 51 was observed (Fig. 4). This is 
consistent with the pattern of a secular resonance 
involving Mercury and Jupiter (plus Venus’ partic¬ 
ipation) (Batygin & Laughlin 2008; Laskar 2008; 
Lithwick & Wu 2011; Bone et al. 2012) and ap¬ 
pears to be the common cause for cm increases in 
nearly all 500-Myr runs, regardless of integrator 
package or coordinates. Note that contributions 
from general relativity (Einstein 1916) as included 
here (Mikkola 1998; Soffel 1989) are important as 
they universally move 51 up (by '^0.43” y“^ at 
present) and away from the (51—55) resonance, 
substantially reducing the probability for large 
Bm increases across all simulations (Laskar 2008; 
Laskar & Gastineau 2009) (Fig. 4). 

More importantly, in the 500-Myr runs body¬ 
centric coordinates showed on average a higher 
tendency towards resonance damping and hence 
larger 51 and smaller ej^^, compared to Jacobi 
coordinates. This tendency is significant and 
roughly halves the probability for cm to increase 
beyond 0.53 in the 500-Myr runs (Fig. 2). The 
cause for the statistical differences originating 


6 



Mercury 



Fig. 4. - Mercury’s maximum eccentricity (e^^) vs. eigen- 

frequency g\ in arcsec y“^ (symbols; g\ from Fourier analysis 
of Mercury’s longitude of perihelion, ■071). m6: mercuryG, HNB: 
HNBody, Bod: bodycentric, Jac: Jacobi coordinates. Dashed 
line: Jupiter’s forcing frequency. For the 5-Gyr runs, and 

g\ values shown were determined using output for and zoi 
over the full 5-Gyr time interval. For the 500-Myr runs, 
and g\ values shown were determined only over the final 250- 
Myr interval of the runs. This improves the representation of 
solutions with long-term decline in below the 0.53-start 
value. Otherwise these solutions would all plot at constant 
=0.53 despite large variations in g\. Contributions from 
general relativity (included in all simulations) universally move 
g\ up (by ~0.43” y~^ at present, arrow) and away from the 
{gi — g^) resonance. 

from different integration coordinates as found 
here is not obvious, neither is clear which (if any) 
of the methods provides accurate probability pre¬ 
dictions over 10®-10®-year timescale (see below). 

5. Mercury’s eccentricity and Lyapunov 
times from different simulations 

Starting from slightly different initial condi¬ 
tions (2.73 m offset in Mercury’s initial radial 
distance), the computed difference in Mercury’s 
eccentricity (Aex) between two solutions grows 
slowly over the first 50 Myr or more when using 
the same program and the same integration co¬ 
ordinates (only bodycentric or only Jacobi coordi¬ 
nates, Fig. 5a-e). This is also true for the same ini¬ 
tial conditions, different programs (mercury6 vs. 
HNBody) and same integration coordinates (both 
bodycentric. Fig. 5f). However, using the same 
initial conditions but different integration coordi¬ 


nates (bodycentric vs. Jacobi coordinates), Ae^ 
between two solutions grows quickly over the first 
10 Myr (Fig. 5g-j). This feature was observed 
irrespective of whether the same or different pro¬ 
grams were used, whether general relativity con¬ 
tributions were included or not, and whether the 
initial eccentricity was low or high (initial cm = 
0.21/0.53, see above). 

During the slow initial rise (Fig. 5a-f), Ae^vj is 
usually dominated by polynomial growth and may 
increase linearly with time in a log-log plot (Laskar 
1989; Varadi et al. 2003). The subsequent rapid 
rise after > 50 Myr until Aex ~ cm has been 
attributed to the chaotic nature of the physical 
system (Laskar 1989; Varadi et al. 2003; Batygin 
& Laughlin 2008). However, it is not clear why 
the system’s chaotic behavior starts dominating 
CjVi’s evolution after e.g. ^60 Myr in mercuryS 
with bodycentric coordinates (Fig. 5a) but only af¬ 
ter ~120 Myr in HNBody with Jacobi coordinates 
(Fig. 5c). The initial Aej^it = 0) in both the 
mercuryS and HNBody simulations was ~3xl0“^^. 
The time evolution of the absolute Acm max¬ 
ima of the two curves (Fig. 5a,c) may be fit to 
simple functions assuming linear and exponential 
growth of the initial Ae^ (Fig- 6). The exponen¬ 
tial function fits the rapid growth phase well in 
the mercuryS simulations (Body-Body) with an 
estimated Lyapunov time r ~ 3.2 Myr. How¬ 
ever, the exponential function is a poor fit to the 
HNBody simulations (Jacobi-Jacobi) with an esti¬ 
mated Lyapunov time r ~ 6 Myr (Fig. 6). 

Nevertheless, these estimated Lyapunov times 
are in good agreement with estimates of Lya¬ 
punov exponents from phase space separation of 
two nearby orbits over time (Fig. 7). Note that 
strictly, Lyapunov exponents are derived from the 
solution of the variational equations, rather than 
from the evolution of two system trajectories (e.g. 
Holman & Murray 1996; Tancredi et al. 2001; Mor- 
bidelli 2002). However, the emphasis here is on 
the relative difference between simulations with, 
say, different integration coordinates, while ap¬ 
plying the same method for estimating Lyapunov 
exponents. Importantly, absolute estimates ob¬ 
tained here are consistent with results from previ¬ 
ous studies (Laskar 1989; Varadi et al. 2003; Baty¬ 
gin & Laughlin 2008). 


7 









m6-m6 LowEcc 6d-6d Body-Body 



0 50 Time (Myr) 100 150 

HNB-HNB LowEcc 4d-4d Jacob-Jacob 



0 50 Time (Myr) 100 150 

HNB-HNB HighEcc 0.25d-0.25d Body-Body 



0 50 Time (Myr) 100 150 

m6-HNB LowEcc 6d-4d Body-Jacob 



m6-m6 HighEcc 2d-2d Body-Body 





0 50 Time (Myr) 100 150 


0 

-2 
-4 
-6 
-8 
-10 

0 50 Time (Myr) 100 150 

HNB-HNB LowEcc 4d-4d Body-Jacob NoPN 


m6-HNB HighEcc 2d-2d Body-Jacob 



h 



Fig. 5. - Computed difference in Mercury’s eccentricity (Ae^) between two runs each per panel for different integrations and 

options. Blue curves: both runs with mercuryG (m6), green curves: both runs with HNBody (HNB), light red curves: one run with 
mercury6, the other with HNBody. Body: bodycentric, Jacob: Jacobi coordinates, md—nd indicates the tiniestep in both runs (in days). 
LowEcc/HighEcc: initial = 0.21/0.53 (see text). NoPN: No contributions from general relativity. Panels a—e show differences in 
between solutions with slightly different initial conditions (2.73 m offset in Mercury’s initial radial distance); panels f—j: identical 
initial conditions. Note rapid Ae>i rise in g—j (all Body-Jacobi). 



























































































































































































































































































































































































log |Ae| Mercury log |Ae| Mercury log |Ae| Mercury 


m6-m6 LowEcc 6d-6d Body-Body 



HNB-HNB LowEcc 4d-4d Jacob-Jacob 



m6-m6 LowEcc 6d-6d Jacob-Jacob NoPN 



Time (Myr) 

6.- Computed difference in Mercury’s eccentricity 

(AeAri) between two runs each per panel with slightly differ¬ 
ent initial conditions (2.73 ni offset in Mercury’s initial radial 
distance) for different integrations and options. Blue curves: 
both runs with mercuryS (m6), green curve: both runs with 
HNBody (HNB), red curves: simple fit functions assuming linear 
and exponential growth of initial Ae^. Body: bodycentric, 
Jacob: Jacobi coordinates, md—nd indicates the timestep in 
both runs (in days). LowEcc: initial = 0.21 (see text). 
NoPN: No contributions from general relativity. Note different 
slopes (estimated Lyapunov time r) of exponential fits in (a) 
and (b) (red dashed lines; linear on logarithmic y-scale). 

The Lyapunov exponents were determined from 
two runs each with mercury6 and HNBody (one 
fiducial, one shadow orbit each) without renormal¬ 
ization and an initial separation of 5x10“^^ AU 
in Mercury’s x-coordinate. Renormalization was 
tested but yielded spurious results as reported be¬ 
fore (Tancredi et al. 2001). If during each time 
interval St, the distance in phase space grows ex¬ 
ponentially, then dj = exp (7 • St), where d^ and 
dj is the initial and final distance. An average jk 
after t = k ■ St may then be computed as: 

A) 

= ’ ( 4 ) 

where St was set to 20 kyr. A log-log plot of 7 ^ 



Fig. 7.- Lyapunov exponents estimated from phase space 

separation of two nearby orbits over time. Blue curve: both 
runs with mercury6 (m6), green curve: both runs with HNBody 
(HNB). Graphs have been truncated shortly after the inflection 
point (see text). Bod: bodycentric, Jac: Jacobi coordinates, 
md—nd indicates timestep in both runs (in days). LowEcc: 
initial ej^ — 0.21. 

vs. time yields a straight line with negative con¬ 
stant slope until an inflection point is reached af¬ 
ter which approaches a non-zero value (Fig. 7). 
The inflection point corresponds to the time in the 
integration where exponential growth starts dom¬ 
inating the evolution of the phase space distance 
between the two orbits. Note that the graphs in 
Fig. 7 have been truncated shortly after the inflec¬ 
tion point to emphasize the plateau in 7 ^ ( 7 ^ keeps 
decreasing subsequently as d is bounded without 
renormalization, not shown). The corresponding 
Lyapunov exponents for the mercuryG and HNBody 
standard runs are and 10“®'^^ (i.e. Lya¬ 

punov times of 2.8 Myr and 5.5 Myr), respectively, 
in good agreement with the Lyapunov times esti¬ 
mated from differences in Mercury’s eccentricity 
evolution (Fig. 6 ). Note that based on different 
approaches it may be difficult to measure Lya¬ 
punov times with an absolute accuracy much bet¬ 
ter than a factor of two (Murray & Holman 1999). 
However, based on the current approach (simu¬ 
lations with identical initial conditions, etc.), it is 
not difficult to measure Lyapunov times with a rel¬ 
ative accuracy better than a factor of two (Figs. 6 
and 7). 

I also computed Ae^vi over time for two simula¬ 
tions with slightly different initial conditions using 
Jacobi coordinates in mercuryB (mixed-variable 
symplectic algorithm) without contributions from 
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general relativity (Fig. 6). In this case, selecting 
parameters for an exponential fit is not obvious 
but it is clear that it takes more than 160 Myr 
for Ae^vf to approach the magnitude of cm- The 
bottom line is that the estimated Lyapunov times 
for Mercury’s orbit from the standard mercuryS 
and HNBody integrations (Figs. 6 and 7) differ at 
least by a factor of two, which is easily measur¬ 
able given the relative accuracy of the current 
approach. However, the chaotic behavior of the 
Solar System is a physical property that cannot 
depend on the numerical algorithm chosen to de¬ 
scribe the system. Yet the results presented here 
suggest that an exponential divergence of trajec¬ 
tories starts dominating the numerical solutions 
after significantly different time intervals, depend¬ 
ing on integration coordinates (Figs. 5, 6, 7). In 
addition, the rapid Ae^vi growth over the first 
10 Myr between two solutions using the same 
initial conditions but different integration coordi¬ 
nates (Fig. 5g-j) suggests a fundamental difference 
between algorithms using bodycentric and Jacobi 
coordinates. 

6. Algorithm tests 

The results presented above raise questions 
about numerical algorithm performance and the 
use of integration coordinates. Why do heliocen¬ 
tric and Jacobi coordinates lead to statistically 
different results? Which algorithm/integration 
coordinates provide more accurate results? This 
section describes several tests aiming at resolv¬ 
ing the dilemma of the dependency of statistical 
results on integration coordinates. 

6.1. 2-Body problem: Position errors 

One approach for testing the accuracy of a nu¬ 
merical algorithm is to study the 2-body prob¬ 
lem, for which an analytical solution exist (strictly, 
Kepler’s equation is solved numerically though). 
While substantially less complex than the general 
n-body problem, in the following the 2-body prob¬ 
lem will illustrate one important difference be¬ 
tween symplectic algorithms with heliocentric and 
Jacobi coordinates. 

The Hamiltonian for the 2-body problem may 
be written as: 

g ^ Irf M! _ _ 

2 mo 2 mi |a;o —a:i| 


where G is the gravitational constant and Pj, m^, 
and Xi are the momenta, masses, and positions of 
the two bodies. The problem can be simplified via 
canonical transformations, using e.g. democratic- 
heliocentric (DH) and Jacobi coordinates, denoted 
here as: 

{x , p) —> {Q , V) : Democratic-Heliocentric 
{x , p) —> {Q , P) : Jacobi . 


In DH coordinates, Qo is the position of the 
center of mass and Qi is the heliocentric posi¬ 
tion of mi. One may use a generating function 
Fm.P) of the new positions and the old mo¬ 
menta (Duncan et al. 1998): 


F3 = 


-Po(So-^Qi) 
-Pi(Qo-^Qi + Qi) , (6) 


where M = mo -I- mi. The relationship between 
old and new variables then is: 


m 

dPo 

dFs 

dPi 

dFs 

dQo 

dF3 

dQi 


xo = {Qo-^Qi) 

xi = (Qo - + 2 i) 

"^0=^0+ Pi 

=Pl - -j^{Po+Pi) ■ 


Note that Pq is the total momentum and Pi is the 
barycentric momentum of mi. After some algebra, 
the Hamiltonian becomes: 


H = 


2 mi 


Gmouii ^ 

) 




2M 


+ 


2mo 


(7) 


The first two terms in parentheses represent the 
Kepler Hamiltonian. The third term represents 
the motion of the center of mass, which moves 
as a free particle and can be ignored. The last 
term needs to be integrated separately in sym¬ 
plectic algorithms with DH coordinates. This 
term is often denoted as FJsun because —Pi = 
Pq — {mo/M){pQ +Pi) is the barycentric momen¬ 
tum of the Sun. 

For Jacobi coordinates, one may use a generat¬ 
ing function 7^2 (®, P) of the old positions and the 
new momenta: 


F 2 = PoimoXo + miXi)/M + Pi{xi - Xq) , ( 8 ) 
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Fig. 8.- Difference between numerically and analytically 

determined position vector of Mercury (|<5r^| = ~ I) 

from 2-body integrations with HNBody of the Sun and Mercury. 
Graphs labeled “Jacobi” and “Body” show |(5r>/t| between the 
2nd order symplectic integrator using Jacobi and bodycentric 
coordinates (both 4-day timestep) and the analytical solution. 
BS: I^T^I between non-symplectic Bulirsch-Stoer algorithm 
(relative accuracy set to 10“^^) and analytical solution. 


where {rriQXQ + miXi)/M = Xs is the center of 
mass. The relationship between old and new vari¬ 
ables then is: 


dF2 

dPo 

dF2 

dPi 

dxo 

dF2 

dxi 


= Qo 

{xi - xq) = Qi 

-Pi + Pomo/M = Pq 
Pi + Pquii/M = Pi 


The Hamiltonian becomes: 


H = 



GrriQmi ^ 

) 



(9) 


where jj, = momi/M is the reduced mass. Again, 
the first two terms represent the Kepler Hamilto¬ 
nian and the third term (center of mass) can be 
ignored. However, comparing the Hamiltonian in 
Jacobi coordinates (Eq. (9)) and in DH coordi¬ 
nates (Eq. (7)) shows that in Jacobi coordinates 
the relevant Kepler mass is p (rather than mi) and 
JJsun is absent. 

Thus, even for the simple 2-body problem, there 
is a difference between symplectic algorithms with 


heliocentric and Jacobi coordinates. Because the 
Hamiltonian is purely Keplerian in Jacobi coor¬ 
dinates (no solar correction), symplectic integra¬ 
tion of the 2-body problem in Jacobi coordinates 
should be more accurate than in DH coordinates. 
Indeed, 10-Myr integrations with HNBody of the 
Sun and Mercury showed that \SrM \ (difference in 
Mercury’s numerical and analytical position vec¬ 
tor) was substantially smaller in integrations with 
Jacobi than DH coordinates over the first 1 Myr or 
so (Fig. 8). On a timescale of 10“^ y, the numeri¬ 
cal position error is close to machine precision for 
Jacobi coordinates but ~10®-times larger for DH 
coordinates (both 4-day timestep). Subsequently, 
the position error grows approximately quadratic 
and linear in time for Jacobi and DH coordinates, 
respectively. Additional runs with HNBody using 
a non-symplectic Bulirsch-Stoer algorithm (rela¬ 
tive accuracy set to 10“^^) yields errors in |rx| 
closer to those of the symplectic integrator with 
Jacobi coordinates (Fig. 8). Note, however, that 
on timescales of 10®-10^ y errors in |r^| grow to 
similarly large values in all three integrations. 

Given that increases in Mercury’s eccentricity 
are critical for the potential destabilization of the 
Solar System, accurate numerical integration of 
its orbit is key. Because Mercury is the inner¬ 
most planet, one might argue that Jacobi coor¬ 
dinates are better suited than DH coordinates for 
this task (Duncan et al. 1998) and that integrating 
Solar System orbits (specifically Mercury’s orbit) 
using DH coordinates simply lacks the necessary 
accuracy. In terms of errors in Mercury’s posi¬ 
tion vector, the present 2-body integrations sup¬ 
port this notion, but only on timescales shorter 
than ~10® y. 

6.2. Bulirsch-Stoer algorithm: Eccentric¬ 
ity errors 

As described above, I also tested whether the 
statistically different in the 500-Myr runs 

(high initial cm ) are related to typically larger er¬ 
rors associated with bodycentric vs. Jacobi coor¬ 
dinates at the same step size. I repeated the 500- 
Myr runs with the bodycentric HNBody setup but 
an eight-fold smaller timestep (0.25 d, ~100-fold 
smaller \AE/E\). The basic result remained the 
same. At high initial cm, the setup using body¬ 
centric coordinates leads to significantly smaller 
mean than the setup using Jacobi coordinates 
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(Fig. 2). This result suggests that the effect of 
different integration coordinates on the statistics 
of Mercury’s eccentricity evolution are not caused 
by errors arising from a too large timestep in the 
bodycentric setup. 

Further tests can be performed by comparing 
the results of the symplectic integrators to results 
of non-symplectic integrators, e.g. the Bulirsch- 
Stoer (BS) algorithm. Note that such tests are 
usually run over shorter time intervals as in most 
cases non-symplectic integrations are computa¬ 
tionally much more expensive than symplectic in¬ 
tegrations. Comparison between HNBody’s BS and 
the 2 nd order symplectic integration of the eight 
planets and Pluto starting at initial = 0.53 
showed a moderate rise in Ae^n over 10 Myr when 
using Jacobi coordinates and a 2-day timestep 
(max|Ae 7 n| ^ 10“^, Fig. 9a). This was not the 
case for bodycentric coordinates/ 2 -day timestep 
vs. BS, where max\AeM\ grew rapidly to ^10“^ 
over 10 Myr (Fig. 9b). However, a moderate rise 
was found again in Ae^n over 10 Myr when using 
bodycentric coordinates and a 0.25-day timestep 
(max|Ae 7 vi| 10“"^, Fig. 9c). Thus, in seem¬ 
ing contradiction to the results of the 500-Myr 
runs (see above), the comparison to BS suggests 
that differences in cm may be rectified provided 
that the timestep in the bodycentric setup is “suf¬ 
ficiently small”. If so, then what is sufficiently 
small? 

For the symplectic 2-body integrations dis¬ 
cussed above, a timestep smaller than ^0.05 days 
would achieve roughly the same accuracy with 
bodycentric coordinates as with Jacobi coordi¬ 
nates and a 4-day timestep. If the same timestep 
ratio were to apply to full Solar System integra¬ 
tions, then the timestep would have to be reduced 
by a factor of 80 when using bodycentric coordi¬ 
nates instead of Jacobi coordinates. For example, 
the 500-Myr runs with 2-day timestep ran over 
roughly 2.2 days wall-clock time. Thus, repeat¬ 
ing those runs with an 80-times reduced timestep 
would take about 180 days (2.6 years for the 5-Gyr 
runs). Not only are such runs currently impracti¬ 
cal, at small timesteps one also needs to consider 
accumulation of numerical errors, which typically 
scale with the number of steps. Finally, could the 
accuracy of symplectic integrations be tested by 
comparison to long-term integrations using non- 
symplectic algorithms such as Bulirsch-Stoer? As 


HNB-HNB HighEcc 2d-1 e-12 Jacob-BS 




HNB-HNB HighEcc 0.25d-1e-12 Body-BS 

-2 1 -^- 1 -^^^^^^^— 



Time (Myr) 


Fig. 9.— Computed difference in Mercury’s eccentricity 
(Ae^) between two runs each per panel for different integra¬ 
tions and options. All runs with HNBody (HNB) including the 
eight planets and Pluto. Body: bodycentric, Jacob: Jacobi co¬ 
ordinates, BS: Bulirsch-Stoer. HighEcc: initial — 0.53 (see 
text), (a) 2nd order symplectic algorithm with 2-day timestep 
and Jacobi coordinates vs. non-symplectic Bulirsch-Stoer al¬ 
gorithm (relative accuracy set to (b) Same as (a) 

but bodycentric coordinates, (c) Same as (b) but 0.25-day 
timestep. 


mentioned above, the Bulirsch-Stoer algorithm is 
currently computationally too expensive for long¬ 
term integrations. In addition, non-symplectic al¬ 
gorithms usually suffer from substantial long-term 
drifts in total energy and angular momentum. 

6.3. Sun + 5 planets: Statistics at high cm 

While long-term integrations of the full equa¬ 
tions of motion of the complete Solar System at 
very small timesteps appear impractical at this 
stage, insight into algorithm performance may be 
gained from test integrations of planetary systems 
of somewhat reduced complexity. It turned out 
that the ( 51 — 55 ) resonance pattern at high initial 
Cm can be reproduced with just a 6 -body setup 
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(Sun + 5 planets: Mercury, Venus, Earth, Jupiter, 
and Saturn). When contributions from general rel¬ 
ativity are ignored, this setup frequently leads to 
rapid cm increases after ~10 Myr. The estimated 
Lyapunov time for Mercury’s orbit in this system 
is ~0.6 Myr. Mercury’s orbit is unstable and in¬ 
termittently switches between resonant and non¬ 
resonant phases (Fig. 10). The resonant phase is 
associated with the (gi— 55 ) resonance and typi¬ 
cally high values in Mercury’s inclination (ca. 10°- 
20°). The non-resonant phase is typically associ¬ 
ated with a drop in Mercury’s inclination and ec¬ 
centricity. After about 8-10 Myr, Acm between 
two nearby orbits reaches the magnitude of cm it¬ 
self, which subsequently increases beyond 0.8 be¬ 
tween 10-12 Myr in some simulation but not in 
others, depending (sensitively) on initial condi¬ 
tions. Hence, the system may provide some useful 
algorithm tests and allow statistical analyses over 
only a 12-Myr interval — an interval short enough 
for integrations with a very small timestep. 

I have integrated 40 solutions each with four dif¬ 
ferent numerical setups over 12 Myr using HNBody 
(Fig. 11): (a) symplectic algorithm with 2-day 
timestep and Jacobi coordinates, (b) symplectic 


Fig. 10.- Mercury’s eccentricity and inclination from 12- 

Myr runs (Sun + 5 planets) at high initial ej^ without con¬ 
tributions from general relativity (NoPN). Jac: Jacobi, Bod: 
bodycentric, coordinates, At: step size in days, BS: Bulirsch- 
Stoer (relative accuracy set to (a) One example from 

each of three setups (see legend, total ^ of runs for each 
setup = 40). (b) Mercury’s inclination for the BS solution 

corresponding to (a). Vertical bars highlight several eccen¬ 
tricity maxima and high inclination values associated with the 
{gi—g^) resonance (BS solution). Arrows indicate a few exam¬ 
ples of eccentricity minima and low inclination values associ¬ 
ated with non-resonant phases. 

algorithm with 0.025-day timestep and bodycen¬ 
tric coordinates, (c) symplectic algorithm with 
0.025-day timestep and Jacobi coordinates, (d) 
Bulirsch-Stoer (BS) algorithm (relative accuracy 
set to 10“^®). All simulations ignore contributions 
from general relativity; the symplectic integrators 
are 2nd order schemes. The setups (b) and (c) 
used an exact timestep of 0.025390625 d (finite bi¬ 
nary representation to minimize round-off errors), 
which is about 80 times smaller than the 2-day 
timestep of setup (a) (see discussion Section 6.2). 
The four sets of integrations {N = 40 each) started 
from the same set of initial conditions, where Mer¬ 
cury’s initial radial distance was offset by 10 m 
between every two adjacent orbits (largest overall 
offset: 39 X 10 m = 390 m). 
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HNB, 5 planets, 12 Myr At = 2d Jac NoPN 


HNB, 5 planets, 12 Myr At = 0.025d Bod NoPN 




Fig. 11.- Mercury’s maximum eccentricity (e^^) achieved 

during final 2 Myr of 12-Myr runs (Sun + 5 planets, see text). 
Jac: Jacobi, Bod: bodycentric, coordinates, At: step size 
in days, BS: Bulirsch-Stoer (relative accuracy set to 
NoPN: No contributions from general relativity. i?o .683 indi¬ 
cates the number of solutions in which grows beyond 0.683 
(mean of all 160 runs, dashed lines). 

For all simulations, Mercury’s maximum eccen¬ 
tricity achieved during the final 2 Myr of 

the 12-Myr runs was recorded, providing a met¬ 
ric for either an increase or a decline in to¬ 
wards the end of the integration (Fig. 10). Of the 
four numerical setups, only the symplectic inte¬ 
grator with At = 0.025 d/Jacobi coordinates and 
the Bulirsch-Stoer algorithm yield similar statisti¬ 
cal results for (Fig. 11). In comparison, the 
symplectic integrator with At = 2 d/Jacobi co¬ 
ordinates shows a reduced tendency for large cm 
increases and has a lower mean value. Strik¬ 
ingly, the symplectic integrator with bodycentric 
coordinates predicts a decline in to values 
below 0.65 between 10-12 Myr in all simulations 
(0.65 is roughly equal to the value from 0- 
10 Myr, see Fig. 10). 

For this particular system, the statistical re¬ 
sults shown in Fig. 10 have important implica¬ 
tions. Based on the agreement between BS and the 
symplectic integrations with At = 0.025 d/Jacobi 
coordinates, it is likely that these two setups pro¬ 
vide the most accurate results. The first impli¬ 
cation then is that At = 2 d is too large for a 
2nd order symplectic integrator with Jacobi coor¬ 


dinates at high Cm- For instance, the calculated 
odds for to increase beyond 0.683 (mean of 
all 160 runs) are 48% at At = 0.025 d but only 
35% at At = 2 d. Second, the symplectic integra¬ 
tor with bodycentric coordinates yields incorrect 
results for this system even at a very small time 
step of 0.025 d = 36 min! The tendency to un¬ 
derestimate increases in Mercury’s eccentricity at 
high Cm , which was also observed in the 500-Myr 
integrations of the full Solar System (Fig. 5), is 
therefore unlikely related to errors associated with 
the size of the timestep. Fundamental differences 
in the implementation of different integration co¬ 
ordinates in symplectic schemes are more likely to 
be the cause (Section 6.1). 

The results of the 12-Myr runs also help resolv¬ 
ing the apparent contradiction between the 500- 
Myr statistics and the Cm comparison to BS men¬ 
tioned in Section 6.2. On the one hand, the 500- 
Myr results suggest that the different statis¬ 
tics are largely independent of the step size in the 
bodycentric setup (Fig. 2). On the other hand, 
the comparison to BS suggests that differences in 
Cm can be minimized when the step size in the 
bodycentric setup is reduced (Fig. 9). In fact, 
for the 12-Myr runs both statements are correct. 
At a small timestep of 0.025 d, Ae^ (BS minus 
Jacobi setup) and Acm (BS minus bodycentric 
setup) are virtually identical during the first 4 Myr 
(not shown). In this interval, Acm is dominated 
by polynomial growth. During the final 2 Myr, 
however, all runs of the bodycentric setup show 
a decline in (contrary to the BS and Jacobi 
setup). In this interval, Acm is dominated by ex¬ 
ponential growth. Importantly, the Hnal interval 
determines the statistics of Mercury’s ultimate ec¬ 
centricity evolution. Thus, a smaller step size does 
improve the bodycentric setup’s accuracy during 
polynomial growth of Acm but has little effect 
during exponential growth of Ae^. 

Note that the implications outlined above 
strictly only apply to the system of reduced com¬ 
plexity studied in this section. That is, a system 
comprising the Sun and just five planets with spe¬ 
cific initial conditions (including high initial e^vi), 
integrated over only 12 Myr and ignoring Post- 
Newtonian corrections. 


14 






















































7. Conclusions 

Reliable predictions of the Solar System’s dy¬ 
namic stability over billions of years not only re¬ 
quire integrators that are fast and accurate but 
also produce robust statistical results in ensem¬ 
ble integrations. Ensemble integrations are neces¬ 
sary because of the chaotic behavior of the system. 
Currently, symplectic integrators are probably the 
best, if not the only, choice in terms of speed and 
accuracy. However, the present results show that 
tackling statistics is trickier as, for instance, the 
predicted probability for a large increase in Mer¬ 
cury’s eccentricity (e^) depends on the choice of 
integration coordinates in symplectic algorithms. 
Several tests performed here suggest that using Ja¬ 
cobi coordinates in symplectic integrations of the 
Solar System is more reliable than using bodycen¬ 
tric coordinates when Mercury’s orbital eccentric¬ 
ity is high. However, these tests reveal little about 
the statistics of long-term integrations when Mer¬ 
cury’s initial orbital eccentricity is low. In fact, 
the influence of Jacobi and bodycentric coordi¬ 
nates on the statistics of Mercury’s eccentricity 
evolution appears to be opposite over 5 Gyr at 
low initial than over 500 Myr at high initial 
Cm- Moreover, even accepting superiority of Ja¬ 
cobi coordinates at high what is the proper 
step size to obtain robust statistical results? If 
applicable, the results of the system of reduced 
complexity (Section 6.3) suggest a timestep that 
could be as small as 0.025 d. Even when starting 
with a larger timestep and reducing the timestep 
during the symplectic integration (which should 
be avoided), such small timesteps would still pose 
a challenge in terms of integration time. Then 
do symplectic integrations with Jacobi coordinates 
and typical timesteps of several days underesti¬ 
mate the odds for disaster (Fig. 11)? Again, if 
applicable, such simulations would underestimate 
the odds for disaster once cm has already reached 
values >0.5 or so. But it is not clear whether the 
odds are predicted correctly to reach those values 
in the first place when starting at low cm- 

It appears that several centuries of analytical 
research, modern state-of-the-art numerical algo¬ 
rithms, and current CPU power has brought us 
closer to answering the question of the Solar Sys¬ 
tem’s long-term stability. However, a definite, ro¬ 
bust answer still seems to be lacking, including an¬ 


swers based on statistical approaches. It is likely 
that the odds for the catastrophic destabilization 
of the inner planets are in the order of a few per¬ 
cent. But what percentage exactly? 
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A. Contributions from general relativity 

Relativistic corrections are critical as they substantially reduce the probability for Mercury’s orbit to 
achieve large eccentricities (Laskar & Gastineau 2009). General relativity (GR) corrections are available 
in HNBody but not in mercuryS. Post-Newtonian (PN) corrections for symplectic integration (Mikkola 
1998; Soffel 1989) were therefore implemented before using mercuryG. For more information on symplectic 
algorithms, see Wisdom & Holman (1991); Saha & Tremaine (1992); Chambers (1999). Note that the PN 
code was fully embedded in the symplectic hybrid integrator mdt_hy(), and not just added as an auxiliary 
force term in mfo_user(). For every planet, one needs to solve an equation of the form (Mikkola 1998): 

a = F = F-k/o(r,t)-k^ /i(r,t!) , (Al) 

where F is the two-body acceleration, is the Hamiltonian perturbation, c is the speed of light, and 

fi{r,v) is the expression for the first Post-Newtonian contribution (Soffel 1989). Combining the implicit 
midpoint method and generalized leapfrog, and denoting the timestep as h, one obtains for the velocity 
jumps (Mikkola 1998): 


5v = hf^ir, t) + ^ (r, v + ^6v) , 


which must be solved for 6v. The first PN acceleration term may be written as (Soffel 1989): 


&■ &■ 


-T* + 4- 

^3 ^ ^4 


r -I- 4-^— 5 -^ V 


Inserting into Eq. (A2) yields: 

bv = /i/o(r,t) 


+h 


GM 


{v+^Sv)'^ GM J’" ■ (^ + , 1C N 

-1 —- r -t 4^- r + M ^^ {v + iJv) 


(A2) 


(A3) 


(A4) 

(A5) 


which was solved iteratively for 6v. Because the relativistic term is small, convergence is rapid (usually 2 
iterations), and [i5t;(i = 2) — 5v[j\/5v[) < 5x10“^® for Mercury. The computational overhead for including GR 
contributions in mercuryG as described above was about 15-20% (wall-clock time). Over the 21st century. 
Mercury’s average perihelion precession (only due to GR) was 0.42976” computed with HNBody and 
0.42978” y“^ computed with mercuryG and the above GR implementation. In terms of Mercury’s long¬ 
term eccentricity {cm) evolution, the difference in cm between mercuryG and HNBody runs (both with GR 
correction. At = 6 days, and bodycentric coordinates) was '^lO”^ over the hrst 40 Myr. 
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Table 2: Initial conditions of the eight planets and 
Pluto“ for 5-Gyr runs from DE431. Heliocentric 
positions r (AU) and velocities v (AU d“^). 



X 

y 

Z 

r 

-1.40712354144735680E-01 

Mercury 

-4.43906230277241465E-01 

-2.33474338281349329E-02 

V 

+2.11691765462179472E-02 

-7.09701275933066148E-03 

-2.52278032052283448E-03 

r 

-7.18629835259113170E-01 

Venus 

-2.25188858612526514E-02 

+4.11716131772919824E-02 

V 

+5.13955712094533914E-04 

-2.03061283748202266E-02 

-3.07198741951420558E-04 

r 

-1.68563248623229384E-01 

Earth + Moon 

+9.68761420122898564E-01 

-1.15183154209270563E-06 

V 

-1.72299715055074729E-02 

-3.01349780674632205E-03 

+2.41254068070491868E-08 

r 

+1.39036162161402177E+00 

Mars 

-2.09984400533893799E-02 

-3.46177919349353047E-02 

V 

+7.47813544105227729E-04 

+1.51863004086334515E-02 

+2.99756038504512547E-04 

r 

+4.00345668418424960E+00 

Jupiter 

+2.93535844833712467E+00 

-1.01823217020834328E-01 

V 

-4.56348056882991196E-03 

+6.44675255807273997E-03 

+7.54565159392195741E-05 

r 

+6.40855153734800886E+00 

Saturn 

+6.56804703677062207E+00 

-3.69127809402511886E-01 

V 

-4.29112154163879215E-03 

+3.89157880254167561E-03 

+1.02876894772680478E-04 

r 

+1.44305195077618524E+01 

Uranus 

-1.37356563056406209E+01 

-2.38128487167790809E-01 

V 

+2.67837949019966498E-03 

+2.67244291355153403E-03 

-2.47764637737944378E-05 

r 

+1.68107582839480649E+01 

Neptune 

-2.49926499733276124E+01 

+1.27271208982211476E-01 

V 

+2.57936917068014599E-03 

+1.77676956230748452E-03 

-9.59089132565213410E-05 

r 

-9.87686582399026491E+00 

Pluto 

-2.79580297772433077E+01 

+5.85080284687055574E+00 

V 

+3.02870206449818878E-03 

-1.53793257901232473E-03 

-7.12171623386267461E-04 


“Masses (Mercury to Pluto in solar masses): 1.66013679527193035E-07, 2.44783833966454472E-06, 3.04043264626852573E- 
06, 3.22715144505387430E-07, 9.54791938424322164E-04, 2.85885980666102893E-04, 4.36625166899970042E-05, 

5.15138902053549668E-05, 7.40740740740740710E-09. 
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